#!/bin/gawk -f

# Script to convert a Biosym mdf file into an Aqua AtomLIB file,
# replacing residue numbers by the wildcard symbol, and 
# keeping only one of every residue.
#
# COMMENTS: AQUA nomenclature based on:
# COMMENTS: Recommendations for presentation of NMR structures of proteins 
# COMMENTS:  and nucleic acids
# COMMENTS: Publication (11/29/96) of Task Group to advice the IUPAC. 
# COMMENTS: Not accepted as of Wed Jul  3 13:40:19 WDT 1996
# COMMENTS: Author: JF Doreleijers
#
# arguments: 
#    type (biosym|akke|disgeo|disman|diana|xplor|pdb|pdbn|type)
#
# mdf format:
# THR_82:CA           C  ca  pepN 0  0   0.1200 0 0 4  1.0000  0.0000 HA CB C N 
#
# aqua format:
# def THR * CA  CA
# def LEU * HB2 HB3
# def PHE * HRX QR
# The first name is in the format determined by 'type', the second is in
# standard internal format, see abovein.

BEGIN {
    if ( type == "" ) type = "biosym"
    type = tolower( type )

    "/bin/date" | getline date
    close "/bin/date"
    "printenv AQUAVERS" | getline vers
    close "printenv AQUAVERS"
    tempfil = "mdf2.tmp"
    tempfil2 = "mdf2.tmp2"
    inmol = 0
    oldres = ""
    oldresnum = 0
    count = 0
    pscount = 0
    rescount = 0
}
function setRS( resnam, atnam ) {
    iupac = atnam
# Code only those which change from Insight's mdf to IUPAC.
#
# Assumes only standard residues are present, look into the code to check
#  wether this will hold for your fancy AA's. OK, ACE of xplor included
# Simple to do it's ONLY CONVERTING 1/2 NOMENCLATURE TO R/S.
#
# NOTE also the heavy atoms of VAL & LEU should;-) could be done
#
# alfa's (only gly not ACE of xplor)
     if (! ( match( resnam, /ACE/ ) ) ) {
      if ( match( atnam, /HA1/ ) ) return ( "HAR" )
      if ( match( atnam, /HA2/ ) ) return ( "HAS" )
     }
# beta's
     if ( match( atnam, /HB/ ) ) { #OK I should have matched the others..
      if ( match(resnam, /PRO/ ) || match(resnam, /GLN/ ) || 
           match(resnam, /GLU/ ) || match(resnam, /LEU/ ) ||
	   match(resnam, /ARG/ ) ||
	   match(resnam, /LYS/ ) || match(resnam, /PHE/ ) ||
           match(resnam, /TRP/ ) || match(resnam, /TYR/ ) ) {
	if ( match( atnam, /HB1/ ) ) return ( "HBR" )
	if ( match( atnam, /HB2/ ) ) return ( "HBS" )
      } else {
        if ( ! (match(resnam, /ALA/ ) ) ) {
	    if ( match( atnam, /HB1/ ) ) return ( "HBS" )
	    if ( match( atnam, /HB2/ ) ) return ( "HBR" )
	}
      }
     } 
# gamma's
     if ( match( atnam, /HG/ ) ) {
      if ( match( atnam, /^HG[12]$/ ) ) {
       if ( match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /GLN/ ) || match( resnam, /GLU/ ) || 
            match( resnam, /CYS/ ) || match( resnam, /MET/ ) ) { 
	if ( atnam == "HG1" ) return ( "HGS" ) 
	if ( atnam == "HG2" ) return ( "HGR" )
       } else { # LYS
        if (match( resnam, /THR/)) return ( "HG1" )
	if ( atnam == "HG1" ) return ( "HGR" ) 
	if ( atnam == "HG2" ) return ( "HGS" )
       }
      }
      if (match(resnam, /ILE/)) {
	if ( atnam == "HG11" ) return ( "HG1R" ) 
	if ( atnam == "HG12" ) return ( "HG1S" )
      }
     }
# delta's
     if ( match( atnam, /HD/ ) ) {
      if (  match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /LYS/ ) )                               {
	if ( match( atnam, /HD1/ ) ) return ( "HDS" )
	if ( match( atnam, /HD2/ ) ) return ( "HDR" )
      } # No exceptions
     }
# epsilon's
     if ( match( atnam, /HE/ ) ) {
      if ( match( resnam, /LYS/ ) )    # Will never occur I guess 
      {
	if ( match( atnam, /HE1/ ) ) return ( "HES" )
	if ( match( atnam, /HE2/ ) ) return ( "HER" )
      }
     }
# if not then just return the proton.
     return ( iupac )
}


function setiupac( resnam, atnam ) {
    iupac = atnam
# Code only those which change from Insight's mdf to IUPAC.
#
# Assumes only standard residues are present, look into the code to check
#  wether this will hold for your fancy AA's. OK, ACE of xplor included
# Simple to do it's ONLY CONVERTING 1/2 NOMENCLATURE TO 2/3.
#
# Protons only:
    if ( match( atnam, /^H/) && ( ! (match( atnam, /X$/ )) ) ) {
#
# alfa's (only GLY not ACE of xplor)
     if (! ( match( resnam, /ACE/ ) ) ) {
      if ( match( atnam, /HA1/ ) ) return ( "HA2" )
      if ( match( atnam, /HA2/ ) ) return ( "HA3" )
     }
# beta's
     if ( match( atnam, /HB/ ) ) {
      if (!( ( match( resnam, /GLY/ ) || match( resnam, /ALA/ ) || 
               match( resnam, /ILE/ ) || match( resnam, /THR/ )  ) ) ) {
	if ( match( atnam, /HB1/ ) ) return ( "HB2" )
	if ( match( atnam, /HB2/ ) ) return ( "HB3" )
      }
     } 
# gamma's
     if ( match( atnam, /HG/ ) ) {
      if ( match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /ASN/ ) || match( resnam, /ASP/ ) || 
            match( resnam, /CYS/ ) || match( resnam, /GLN/ ) || 
            match( resnam, /GLU/ ) || match( resnam, /ILE/ ) || 
            match( resnam, /LYS/ ) || match( resnam, /MET/ ) ) {
        
	if ( atnam == "HG1" ) return ( "HG2" ) 
	if ( atnam == "HG2" ) return ( "HG3" )
	if ( atnam == "HG11" ) return ( "HG12" ) # For Ile case
	if ( atnam == "HG12" ) return ( "HG13" )
      }
     }
# delta's
     if ( match( atnam, /HD/ ) ) {
      if ( match( resnam, /PRO/ ) || match( resnam, /ARG/ ) || 
            match( resnam, /GLN/ ) || match( resnam, /GLU/ ) ||
            match( resnam, /LYS/ ) )                               {
	if ( match( atnam, /HD1/ ) ) return ( "HD2" )
	if ( match( atnam, /HD2/ ) ) return ( "HD3" )
      }
     }
# epsilon's
     if ( match( atnam, /HE/ ) ) {
      if ( match( resnam, /LYS/ ) )                                 {
	if ( match( atnam, /HE1/ ) ) return ( "HE2" )
	if ( match( atnam, /HE2/ ) ) return ( "HE3" )
      }
     }
# if not then just return the proton.
     return ( iupac )
    }

#
# Pseudo's only. Always H?X -> Q? in case of methylenes
#                    so just do the methyls and ILE QG and a few others
    if ( match( atnam, /^H/ ) && match ( atnam, /X$/ ) ) {
# alfa's (no exceptions)
# beta's 
     if ( match( atnam, /HB/ ) && match( resnam, /ALA/ ) ) return ( "MB" )
# gamma's
     if ( match( atnam, /HG/ ) ) {
      if ( match( resnam, /VAL/ ) ) {
	if ( match( atnam, /HG1/ ) ) return ( "MG1" )
	if ( match( atnam, /HG2/ ) ) { return ( "MG2" ) 
        } else                         return ( "QG" )
      }
      if ( match( resnam, /ILE/ ) ) {
	if ( match( atnam, /HG1/ ) ) return ( "QG" )
	if ( match( atnam, /HG2/ ) ) return ( "MG" )
      }
      if ( match( resnam, /THR/ ) ) {
	if ( match( atnam, /HG2/ ) ) return ( "MG" )
      }
     }
# delta's
     if ( match( atnam, /HD/ ) ) {
      if ( match( resnam, /LEU/ ) ) {
	if ( match( atnam, /HD1/ ) ) return ( "MD1" )
	if ( match( atnam, /HD2/ ) ) { 
            return ( "MD2" ) 
	} else                         return ( "QD" )
      }
      if ( match( resnam, /ILE/ ) )    return ( "MD" )
      if ( match( resnam, /ASN/ ) )    return ( "QD" )
     }
# epsilon's
     if ( match( atnam, /HE/ ) ) {
      if ( match( resnam, /MET/ ) ) return ( "ME" )
      if ( match( resnam, /GLN/ ) ) return ( "QE" )
     }
# zeta's
     if ( match( atnam, /HHX/ ) && match( resnam, /ARG/ ) ) return ( "QZ" )
# what remains: Always H?X -> Q? this covers methylenes, arg amino, rings
     dummy = match( atnam, /X$/ ) # Get the RSTART value
     at = substr( atnam, 2, RSTART-2 )
     iupac = sprintf( "Q%s", at )
     return ( iupac )
    }
# heavy atoms are all the same so just return it.
    return ( iupac )
}

function setpdb( resnam, atnam,    seq, at ) {
    pdb = atnam
    if ( atnam == "HN" ) 
	return ( "H" )
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /[123]$/ ) ) {
	    seq = substr( atnam, RSTART, 1 )
	    at = substr( atnam, 1, RSTART-1 )
	    if ( at != "HB" ) {
		if ( match( resnam, /TYR/ ) || match( resnam, /PHE/ ) || 
		     match( resnam, /HIS/ ) || match( resnam, /TRP/ ) ) 
		    return ( pdb )
		if ( at == "HG" && match( resnam, /THR/ ) )
		    return ( pdb )
	    }
	    pdb = sprintf( "%s%s", seq, at )
	}
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    pdb = sprintf( "Q%s", at )
	}
    }
    if ( pdb == "QG" )
	if ( match( resnam, /VAL/ ) )
	    pdb = "QQG"
    if ( pdb == "QD" )
	if ( match( resnam, /LEU/ ) )
	    pdb = "QQD"
    return ( pdb )
}

function setpdbn ( resnam, atnam,    seq, at ) {

# The same as setpdb but now with the '2/3 flavor'

    pdbn = atnam
    
    if ( atnam == "HN" ) 
	return ( "H" )
	
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /[123]$/ ) ) {
	
	    atnam = setiupac( resnam, atnam )          # The magic trick..
	    match( atnam, /[123]$/ )                   # to get RSTART again
	    
	    seq = substr( atnam, RSTART, 1 )
	    at = substr( atnam, 1, RSTART-1 )
	    if ( at != "HB" ) {
		if ( match( resnam, /TYR/ ) || match( resnam, /PHE/ ) || 
		     match( resnam, /HIS/ ) || match( resnam, /TRP/ ) ) 
		    return ( pdbn )
		if ( at == "HG" && match( resnam, /THR/ ) )
		    return ( pdbn )
	    }
	    
	    pdbn = sprintf( "%s%s", seq, at )
	}
	
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    pdbn = sprintf( "Q%s", at )
	}

    }

    if ( pdbn == "QG" )
	if ( match( resnam, /VAL/ ) )
	    pdbn = "QQG"
    if ( pdbn == "QD" )
	if ( match( resnam, /LEU/ ) )
	    pdbn = "QQD"

    return ( pdbn )
}


function setdiana( resnam, atnam ) {
    diana = setdisgeo( resnam, atnam )
    return ( diana )
}


function setdisman( resnam, atnam,    aro ) {
    disman = setdisgeo( resnam, atnam )
    if ( resnam == "HIS" || resnam == "PHE" || resnam == "TYR" || resnam == "TRP" )
	aro = setnmraro( resnam, atnam )
    else
	aro = ""
    if ( aro != "" )
	return ( aro )
    else
	return ( disman )
}

function setdisgeo( resnam, atnam,    at ) {
    disgeo = atnam
    if ( atnam == "HN" ) 
	return ( "H" )
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    disgeo = sprintf( "Q%s", at )
	}
    }
    if ( disgeo == "QG" )
	if ( match( resnam, /VAL/ ) )
	    disgeo = "QQG"
    if ( disgeo == "QD" )
	if ( match( resnam, /LEU/ ) )
	    disgeo = "QQD"
    return ( disgeo )
}

function setakke( resnam, atnam, mult,   at ) {
    akke = atnam
    if ( match( atnam, /^H/ ) ) {
	if ( match( atnam, /X$/ ) ) {
	    at = substr( atnam, 2, RSTART-2 )
	    if ( mult == 3 )
		akke = sprintf( "M%s", at )
	    else
		akke = sprintf( "QH%s", at )
	}
    }
    if ( match( resnam, /ILE/ ) ) {
	if ( akke == "QHG1" )
	    akke = "QHG"
	if ( akke == "MD1" )
	    akke = "MD"
	if ( akke == "MG2" )
	    akke = "MG"
	if ( akke == "MD1" )
	    akke = "MD"
    }
    if ( match( resnam, /THR/ ) )
	if ( akke == "MG2" )
	    akke = "MG"
    if ( match( resnam, /ASN/ ) )
	if ( akke == "QHD2" )
	    akke = "QHD"
    if ( match( resnam, /GLN/ ) )
	if ( akke == "QHE2" )
	    akke = "QHE"
    if ( match( resnam, /VAL/ ) )
	if ( akke == "MG" )
	    akke = "QMG"
    if ( match( resnam, /LEU/ ) )
	if ( akke == "MD" )
	    akke = "QMD"
    if ( match( resnam, /TYR/ ) || match( resnam, /PHE/ ) ) {
	if ( akke == "QHD" )
	    akke = "HD"
	if ( akke == "QHE" )
	    akke = "HE"
    }
    return ( akke )
}

function setxplor( resnam, atnam ) {
    xplor = atnam
    if ( atnam == "HRX" ) return ( "skip" )
    if ( match( xplor, /^H/ ) )
	sub( /X$/, "#", xplor ) 
#    if ( xplor == "HG#" )
#	if ( match( resnam, /VAL/ ) )
#	    xplor = "HG*"
#    if ( xplor == "HD#" )
#	if ( match( resnam, /LEU/ ) )
#	    xplor = "HD*"
    return ( xplor )
}

function setnmraro( resnam, atnam,    at, atext ) {
# set the nmr numbered name for ring atoms
    at = substr( atnam, 1, 1)
    atext = atnam
    sub( at, "", atext )
    if ( resnam == "HIS" ) {
	if ( atext == "D1" ) return ( at "D1" )
	if ( atext == "E1" ) return ( at "E2" )
	if ( atext == "E2" ) return ( at "E3" )
	if ( atext == "D2" ) return ( at "D4" )
    }
    if ( resnam == "TRP" ) {
	if ( atext == "E1" ) return ( at "E1" )
	if ( atext == "D1" ) return ( at "D2" )
	if ( atext == "G" )  return ( at "G3" )
	if ( atext == "E3" ) return ( at "E4" )
	if ( atext == "Z3" ) return ( at "Z5" )
	if ( atext == "H2" ) return ( at "H6" )
	if ( atext == "Z2" ) return ( at "Z7" )
	if ( atext == "E2" ) return ( at "E8" )
	if ( atext == "D2" ) return ( at "D9" )
    }
    if ( resnam == "TYR" || resnam == "PHE" ) {
	if ( atext == "G" )  return ( at "G1" )
	if ( atext == "D1" ) return ( at "D2" )
	if ( atext == "E1" ) return ( at "E3" )
	if ( atext == "Z" )  return ( at "Z4" )
	if ( atext == "E2" ) return ( at "E5" )
	if ( atext == "D2" ) return ( at "D6" )
        if ( atnam == "HDX" ) return ( "skip" )
        if ( atnam == "HEX" ) return ( "skip" )
    }
    return ( "" )
}

function store( resnam, resnum, atnam ) {
# if the residue number just changed, check whether the
# previous residue has to be stored in the residue list
    if ( resnum != oldresnum ) {
	stored = 0
	for ( i=1; i<=rescount; i++ ) {
	    if ( res[i] == oldres )
		stored = 1
	}
	if ( !stored )
	    res[++rescount] = oldres
	oldresnum = resnum
	oldres = resnam
    }
# has current residue been done already?
    for ( i=1; i<=rescount; i++ ) {
	if ( res[i] == resnam )
	    return
    }
# store residue and atom names
    atr[++count] = resnam
    atn[count] = atnam
}

/^#/ {
    next
}

/@molecule/ {
    mol = $2":"
    inmol = 1
    next
}

/@list pseudoatom/ {
    inmol = 0
    sub( mol, "" )
    sub( /_[0123456789]*:/, " " )
# check whether this pseudo atom occurred before
# and skip some special cases
    for ( i=1; i<=pscount; i++ ) {
	if ( psr[i] == $3 && psn[i] == $4 )
	    next
    }
    if ( $3 == "ARG+" || $3 == "SERN" ) next
    if ( $3 == "ALAn" || $3 == "VALn" ) next
    if ( $3 == "CA" || $3 == "OH" ) next
# store the pseudoatom residue and atom names
    psr[++pscount] = $3
    psn[pscount] = $4
# store the number of constituent atoms (assuming they are on one line)
    getline
    psc[pscount] = NF
}

inmol {
    if ( NF == 0 ) 
	next
    inmol = !match( $1, /^@/ )
    if ( !inmol ) 
	next
    sub( /_/, " " )
    sub( /:/, " " )
# skip some special cases
    if ( $1 == "ARG+" || $1 == "SERN" ) next
    if ( $1 == "ALAn" || $1 == "VALn" ) next
    if ( $1 == "CA" || $1 == "OH" ) next
# store the atom
    store( $1, $2, $3 )
}


END {
# check whether the last residue has to be stored in the residue list
# (cf. function store)
    stored = 0
    for ( i=1; i<=rescount; i++ ) {
	if ( res[i] == oldres )
	    stored = 1
    }
    if ( !stored )
	res[++rescount] = oldres

# print the residue list, sort, and read it again
    for ( i=1; i<=rescount; i++ )
	print res[i] > tempfil
    rescount = 0
    while ( "/bin/sort "tempfil | getline )
	res[++rescount] = $1
    close "/bin/sort "tempfil
    system( "/bin/rm "tempfil )

# output the atoms in the order of the ordered residue list
    total = 0
    for ( i=1; i<=rescount; i++ ) {
	inres = 0
	for ( j=1; j<=count; j++ ) {
	    if ( atr[j] == res[i] ) {
		inres = 1
		intern = setiupac( atr[j], atn[j] )
		if ( type == "type" )
		    printf( "def %-5s * %-5s\n", atr[j], intern ) > tempfil2
		else {
		    if ( type == "pdb" )    extern = setpdb( atr[j], atn[j] )
		    if ( type == "pdbn" )   extern = setpdbn( atr[j], atn[j] )
		    if ( type == "akke" )   extern = setakke( atr[j], atn[j], 1 )
		    if ( type == "disgeo" ) extern = setdisgeo( atr[j], atn[j] )
		    if ( type == "disman" ) extern = setdisman( atr[j], atn[j] )
		    if ( type == "diana" )  extern = setdiana( atr[j], atn[j] )
		    if ( type == "xplor" )  extern = setxplor( atr[j], atn[j] )
		    if ( type == "biosym" ) extern = atn[j]
		    if ( extern != "skip" ) {
			printf( "def %-5s * %-5s %-5s\n", atr[j], extern, intern ) > tempfil2
			total++
		    }
		    else
			continue
		    if ( type == "biosym" && match( atn[j], /^H.*[123]$/ ) ) {
			new = setRS(atr[j], atn[j])
			if (new != atn[j]) {
				printf( "def %-5s * %-5s %-5s\n", atr[j], new, intern ) > tempfil2
				total++
			}
		    }
		    
		    if ( ( type == "pdb") && match( atn[j], /^H.*[123]$/ ) ) {
			printf( "def %-5s * %-5s %-5s\n", atr[j], atn[j], intern ) > tempfil2
			total++
		    }
		    
		    if ( ( type == "pdbn") && match( atn[j], /^H.*[123]$/ ) ) {
			printf( "def %-5s * %-5s %-5s\n", atr[j], intern, intern ) > tempfil2
			total++
		    }
		    
		    if ( ( (type == "pdb") || (type == "pdbn") ) && atn[j] == "HN" ) {
			printf( "def %-5s * %-5s %-5s\n", atr[j], "HN", "HN" ) > tempfil2
			total++
		    }
		}
	    }
	    else
		if ( inres ) 
		     break
	}
	for ( j=1; j<=pscount; j++ ) {
	    if ( psr[j] == res[i] ) {
		intern = setiupac( psr[j], psn[j] )
		if ( type == "type" )
		    printf( "def %-5s * %-5s\n", psr[j], intern ) > tempfil2
		else {
		    if ( type == "pdb" )    extern = setpdb( psr[j], psn[j] )
		    if ( type == "pdbn" )   extern = setpdbn( psr[j], psn[j] )
		    if ( type == "akke" )   extern = setakke( psr[j], psn[j], psc[j] )
		    if ( type == "disgeo" ) extern = setdisgeo( psr[j], psn[j] )
		    if ( type == "disman" ) extern = setdisman( psr[j], psn[j] )
		    if ( type == "diana" )  extern = setdiana( psr[j], psn[j] )
		    if ( type == "xplor" )  extern = setxplor( psr[j], psn[j] )
		    if ( type == "biosym" ) extern = psn[j]
		    if ( extern != "skip" ) {
			printf( "def %-5s * %-5s %-5s\n", psr[j], extern, intern ) > tempfil2
			total++
		    }
		    else
			continue
		    if ( type == "biosym" && match( psn[j], /^H.*X$/ ) ) {
			new = psn[j]
			sub( /X$/, "*", new )
			printf( "def %-5s * %-5s %-5s\n", psr[j], new, intern ) > tempfil2
			total++
			if ( match(psn[j],/HD[12]/) && match(psr[j],/LEU/)) {
			  if (match(intern,/MD1/)) extern = "HDR*"
			  if (match(intern,/MD2/)) extern = "HDS*"
			  printf( "def %-5s * %-5s %-5s\n", psr[j], extern, intern ) > tempfil2
			  total++
			}
			if ( match(psn[j],/HG[12]/) && match(psr[j],/VAL/)) {
			  if (match(intern,/MG1/)) extern = "HGR*"
			  if (match(intern,/MG2/)) extern = "HGS*"
			  printf( "def %-5s * %-5s %-5s\n", psr[j], extern, intern ) > tempfil2
			  total++
			}
		    }
		}
	    }
	}
    }

# output the header and cat the file
    if ( type == "type" ) {
	print "$ Aqua Version", vers, "DATA: Atom type definitions"
	print "$ FILE: AtomTypeLIB CREATED:", date
    }
    else {
	print "$ Aqua Version", vers, "DATA: Atom name definitions " toupper(type)
	print "$ FILE: AtomLIB CREATED:", date
    }
    print "$ REMARK: "
    print "$ REMARK: AQUA names based on:"
    print "$ REMARK: Recommendations for presentation of NMR structures of proteins "
    print "$ REMARK:  and nucleic acids"
    print "$ REMARK: Publication (11/29/96) of Task Group to advice the IUPAC. "
    print "$ REMARK: Author: JF Doreleijers and JAC Rullmann"
    print "$ REMARK: "
    print "$ ATOM_COUNT:", total
    print "$ BEGIN_DATA:"
    system( "/bin/cat "tempfil2 "; /bin/rm "tempfil2 )

}


